Introduction to Finite Element Analysis 
(FEA) or Finite Element Method (FEM) 


Finite Element Analysis (FEA) or Finite 
Element Method (FEM) 


€ The Finite Element Analysis (FEA) is a 
numerical method for solving problems of 
engineering and mathematical physics. 


€ Useful for problems with complicated 
geometries, loadings, and material properties 
where analytical solutions can not be 
obtained. 


The Purpose of FEA 


Analytical Solution 


e Stress analysis for trusses, beams, and other simple 
structures are carried out based on dramatic simplification 
and idealization: 


— mass concentrated at the center of gravity 
— beam simplified as a line segment (same cross-section) 


* Design is based on the calculation results of the idealized 
structure & a large safety factor (1.5-3) given by experience. 


FEA 


* Design geometry is a lot more complex; and the accuracy 

requirement is a lot higher. We need 

— To understand the physical behaviors of a complex 
object (strength, heat transfer capability, fluid flow, etc.) 

— To predict the performance and behavior of the design; 
to calculate the safety margin; and to identify the 
weakness of the design accurately; and 

— To identify the optimal design with confidence 


Brief History 


Grew out of aerospace industry 

è Post-WW II jets, missiles, space flight 
€ Need for light weight structures 

e Required accurate stress analysis 

e Paralleled growth of computers 
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Common FEA Applications 


Mechanical/Aerospace/Civil/Automotive 
Engineering 
Structural/Stress Analysis 
ㆍ Static/Dynamic 
ㆍ Linear/Nonlinear 
Fluid Flow 
Heat Transfer 
Electromagnetic Fields 
Soil Mechanics 
Acoustics 
Biomechanics 


Discretization 


Complex Object Simple Analysis 
(Material discontinuity, 


Complex and arbitrary geometry) 


Simplified : ! : 

Real Mathematical Discretized 

Word | > | (Idealized) | => Model —> | (mesh) 
Physical 


Model 


Model 


Discretizations 


+ Model body by dividing it into an 
equivalent system of many smaller bodies 
or units (finite elements) interconnected at 
points common to two or more elements 
(nodes or nodal points) and/or boundary 
lines and/or surfaces. 


Types of Finite Elements 


1-D (Line) Element 


(Spring, truss, beam, pipe, etc.) 


3-D (Solid) Element 
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(Membrane, plate, shell, etc.) 


L 
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(3-D fields - temperature, displacement, stress, flow velocity) 


Elements & Nodes - Nodal Quantity 


Feature 


+ Obtain a set of algebraic equations to 
solve for unknown (first) nodal quantity 
(displacement). 


+ Secondary quantities (stresses and 
strains) are expressed in terms of nodal 


values of primary quantity 


Object 


Elements 


Nodes 


A. 


| 


Displacement | ——» 
— 


Examples of FEA — 1D (beams) 


Examples of FEA - 2D 


Examples of FEA — 3D 


Advantages 


Irregular Boundaries 

General Loads 

Different Materials 

Boundary Conditions 

Variable Element Size 

Easy Modification 

Dynamics 

Nonlinear Problems (Geometric or Material) 
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The following notes are a summary from “Fundamentals of Finite Element Analysis” by David V. Hutton 


Principles of FEA 


The finite element method (FEM), or finite element analysis 
(FEA), 1s a computational technique used to obtain approximate 
solutions of boundary value problems in engineering. 


Boundary value problems are also called field problems. The field 
1s the domain of interest and most often represents a physical 
structure. 


The field variables are the dependent variables of interest governed 
by the differential equation. 


The boundary conditions are the specified values of the field 
variables (or related variables such as derivatives) on the boundaries 
of the field. 


For simplicity, at this point, we assume a two-dimensional case with a 
single field variable g(x, y) to be determined at every point P(x, y) such 


that a known governing equation (or equations) is satisfied exactly at every 
such point. 


(a) (b) 
(a) A general two-dimensional domain of field variable (x, y). 


(b) A three-node finite element defined in the domain. (c) Additional 
elements showing a partial finite element mesh of the domain. 


(c) 


-A finite element is not a differential element of size dx X dy. 


- A node 1s a specific point in the finite element at which the value of the 
field variable 1s to be explicitly calculated. 


Shape Functions 


The values of the field variable computed at the nodes are used to 
approximate the values at non-nodal points (that 1s, in the element interior) 
by interpolation of the nodal values. For the three-node triangle example, 
the field variable 1s described by the approximate relation 


0( y) = N(x, y) 07 + Ny, y) 02 + N3(%, y) 95 


where 0), p, and p; are the values of the field variable at the nodes, and 
N,, N,, and N; are the interpolation functions, also known as shape 
functions or blending functions. 


In the finite element approach, the nodal values of the field variable are 
treated as unknown constants that are to be determined. The interpolation 
functions are most often polynomial forms of the independent variables, 
derived to satisfy certain required conditions at the nodes. 


The interpolation functions are predetermined, known functions of the 
independent variables; and these functions describe the variation of the 
field variable within the finite element. 


Degrees of Freedom 


Again a two-dimensional case with a single field variable (x, y). The 
triangular element described is said to have 3 degrees of freedom, as three 
nodal values of the field variable are required to describe the field variable 
everywhere in the element (scalar). 


p(x, y) = N(x, y) 07 + Nox, y) e; + Ni, y) O 


In general, the number of degrees of freedom associated with a finite 
element is equal to the product of the number of nodes and the number of 
values of the field variable (and possibly its derivatives) that must be 
computed at each node. 


A GENERAL PROCEDURE FOR 
FINITE ELEMENT ANALYSIS 


Preprocessing 


Define the geometric domain of the problem. 
— Define the element type(s) to be used (Chapter 6). 
— Define the material properties of the elements. 
— Define the geometric properties of the elements (length, area, and the like). 
— Define the element connectivities (mesh the model). 
Define the physical constraints (boundary conditions). Define the loadings. 


Solution 


— computes the unknown values of the primary field variable(s) 


— computed values are then used by back substitution to compute additional, derived variables, such as 
reaction forces, element stresses, and heat flow. 


Postprocessing 


—  Postprocessor software contains sophisticated routines used for sorting, printing, and plotting 
selected results from a finite element solution. 


Stifíness Matrix 


The primary characteristics of a finite element are embodied in the 
element stiffness matrix. For a structural finite element, the 
stiffness matrix contains the geometric and material behavior 
information that indicates the resistance of the element to 
deformation when subjected to loading. Such deformation may 
include axial, bending, shear, and torsional effects. For finite 
elements used in nonstructural analyses, such as fluid flow and heat 
transfer, the term stiffness matrix 1s also used, since the matrix 
represents the resistance of the element to change when subjected 
to external influences. 


LINEAR SPRING AS A FINITE ELEMENT 


A linear elastic spring 19 a mechanical device capable of supporting axial 
loading only, and the elongation or contraction of the spring 1s directly 
proportional to the applied axial load. The constant of proportionality 
between deformation and load is referred to as the spring constant, spring 
rate, or spring stiffness k, and has units of force per unit length. As an 
elastic spring supports axial loading only, we select an element coordinate 
system (also known as a local coordinate system) as an x axis oriented 
along the length of the spring, as shown. 


Force, f 


uy 7 uy k 
1 
fi 1 k 2 h x Deflection, ô = u, — u) 


(a) (b) 


(a) Linear spring element with nodes, nodal displacements, and nodal forces. 
(b) Load-deflection curve. 


Assuming that both the nodal displacements are zero when the spring is 
undeformed, the net spring deformation 1s given by 
Ô= Uy — Uy 
and the resultant axial force in the spring 18 
f= kò = k(u, — uj) 
For equilibrium, 
Sith =9 or 4 =~ hy 
Then, in terms of the applied nodal forces as 
1 = 04. - uj) 
f; = k(u, — uy) 
which can be expressed in matrix form as 


k — —k|Juil Jf 
Au AS CN rams 


where 


k k 
[ke] = a Stiffness matrix for one spring element 


k KE 


system), {u} is the column matrix (vector) of nodal displacements, and { f } is the 
column matrix (vector) of element nodal forces. 


AL xs Uy _| k —k 
| h | Mò | 162 | is ike] = k | 


The equation shows that the element stiffness matrix for the linear spring element 
isa2 X 2 matrix. This corresponds to the fact that the element exhibits two nodal 
displacements (or degrees of freedom) and that the two displacements are not 
independent (that is, the body is continuous and elastic). 


Furthermore, the matrix is symmetric. This is a consequence of the symmetry of 
the forces (equal and opposite to ensure equilibrium). 


Also the matrix is singular and therefore not invertible. That 1s because the 
problem as defined 1s incomplete and does not have a solution: boundary 
conditions are required. 


SYSTEM OF TWO SPRINGS 


These are internal forces 


Writing the equations for each spring in matrix form: 
AA Superscript refers to element 


(2) 
g | po | 
J 3 
To begin assembling the equilibrium equations describing the behavior of the 


system of two springs, the displacement compatibility conditions, which relate 
element displacements to system displacements, are written as: 


T 1 T 1 
| ox NET 
2T NM ~~ — 
N — 
>= | > | 
No An — e 
N — 
——— | — 
— — 
= = = = 
ore were 
—.= — — — 
| | 


(1) (1) (2) (2) 
uy =U, u, = U2 u, = U2 u, = Us 


And 


(1) 
therefore: ki ki JUL] _ J 71 


Here, we use the notation f'/),to represent the force exerted on element j at node i. 


Expand each equation in matrix form: 


Summing member by member: 
(1) 


서 -h O7fU | 
E ki + ka IPIE fo. yo 


0 —Ko ka U3 po 
J3 


Next, we refer to the free-body diagrams of each of the three nodes: 


fy-n fje-f?7-R  fy-Fh 


Final form: 


ko =k OTU F, 
-k kitke -k Ut =! F, (1) 
0 — ko kò U4 F3 


Where the stiffness matrix: 


ki —ki 0 
kis [ 4 ki + ka I 


0 —k» ka 


Note that the system stiffness matrix 1s: 

(1) symmetric, as is the case with all linear systems referred to orthogonal coordinate 
systems; 

(2) singular, since no constraints are applied to prevent rigid body motion of the 
system; 

(3) the system matrix 19 simply a superposition of the individual element stiffness 
matrices with proper assignment of element nodal displacements and associated 
stiffness coefficients to system nodal displacements. 


FEA for multiple (many) elements { F} = [ K ] . {U} 
ht 


"DT di | Array of dispiscements | cne 
(one for each DOF) Matrix of for each DOF) 
zofinesses 
(DOF x DOF) 


{F ! 1s "known" (loads) 
[K ] is “known” (geometry, material properties. ..elements) 

{U} is to be determined (displacements) 

This can be solved mathematically using a matrix inversion method 
(Fi=[K] t) > w= iF) 


(first nodal quantity) 


Once the displacements {U} are known, then strains and stresses can be determined: 


2m = (1-D ...more complicated for 2-D and 3-D strains) 


c=E-e 


c, T 
and FOS = — (second nodal quantities) 
c 


Example with Boundary Conditions 


Consider the two element system as described before where Node 1 is attached to a 
fixed support, yielding the displacement constraint U, = 0, k,=50 Ib/in, k,=75 Ib/in, 
Fy- F= 75 lb for these conditions determine nodal displacements U, and U}. 


anman. e 
1 fP + fe 2 F, fe F, 


(c) (d) (e) 


Substituting the specified values into (1) we have: 


yo Due to boundary condition 


Example with Boundary Conditions 


Because of the constraint of zero displacement at node 1, nodal force F', becomes an 


unknown reaction force. Formally, the first algebraic equation represented in this 
matrix equation becomes: 


-SOU, =F, 


and this is known as a constraint equation, as it represents the equilibrium condition 
of a node at which the displacement is constrained. The second and third equations 


become 
125 -75]fU,] _ | 75 
-75 75 Jiu] (75 
which can be solved to obtain U, = 3 in. and U} = 4 in. Note that the matrix 
equations governing the unknown displacements are obtained by simply striking out 
the first row and column of the 3 X 3 matrix system, since the constrained 


displacement is zero (homogeneous). If the displacement boundary condition is not 


equal to zero (nonhomogeneous) then this is not possible and the matrices need to be 
manipulated differently (partitioning). 


Truss Element 


The spring element is also often used to represent the elastic nature of supports for 
more complicated systems. A more generally applicable, yet similar, element 1s an 
elastic bar subjected to axial forces only. This element, which we simply call a bar or 
truss element, 1s particularly useful in the analysis of both two- and three- 
dimensional frame or truss structures. Formulation of the finite element 
characteristics of an elastic bar element 19 based on the following assumptions: 


1.The bar is geometrically straight. 

2.The material obeys Hooke’s law. 

3.Forces are applied only at the ends of the bar. 

4.The bar supports axial loading only; bending, torsion, and shear are not 
transmitted to the element via the nature of its connections to other elements. 


Truss Element Stiffness Matrix 


Let's obtain an expression for the stiffness matrix K for the beam element. Recall 
from elementary strength of materials that the deflection 0 of an elastic bar of 
length L and uniform cross-sectional area A when subjected to axial load P : 


PL 
AE 
where £ is the modulus of elasticity of the material. Then the equivalent spring 
constant k: 
AE 
k a ————— 
0 L 


Therefore the stiffness matrix for one element is: 


AE| 1 —] 


And the equilibrium equation in matrix form: 


피그 가바 


Truss Element Blending Function 


An elastic bar of length L to which is affixed a uniaxial coordinate system x with its 
origin arbitrarily placed at the left end. This is the element coordinate system or 
reference frame. Denoting axial displacement at any position along the length of the 
bar as u(x), we define nodes | and 2 at each end as shown and introduce the nodal 
displacements: 


u,=u (x=0) and u, = u(x = L). 

Thus, we have the continuous field variable u(x), which is to be expressed 
(approximately) in terms of two nodal variables u, and u,. To accomplish this 
discretization, we assume the existence of interpolation functions N(x) and N,(x) 


(also known as shape or blending functions) such that 


u(x) = N (x)u, + Ny(x)u; 


E 
n m 


l 
x 이에 m u(x) | 
< L — 


2 x 


Truss Element Blending Function 


To determine the interpolation functions, we require that the boundary values of u(x) 
(the nodal displacements) be identically satisfied by the discretization such that: 


u,—u (x=0) and u, = u(x = L). 
lead to the following boundary (nodal) conditions: 


N,(0)=1 N,(0) =0 
N(L)=0 NXL)-1 


As we have two conditions that must be satisfied by each of two one-dimensional 
functions, the simplest forms for the interpolation functions are polynomial forms: 


N(x) = ay + ax 
Nx) = by + bx 


Truss Element Blending Function 


Where the polynomial coefficients are to be determined via satisfaction of the 
boundary (nodal) conditions. Application of conditions yields a, =1,b,=0, 
therefore a, = —(1/L) and b, = x/L. Therefore, the interpolation functions are 


N(x) =1-x/L 
Nx) = x/L 


total interpolation 


— 


Therefore the final expression of the blending function is: 


u(x) = (1 — x/L)u, + (x/L)u, Uy 


u 
1 
And in matrix form: 


u, contribution 


u(x) = [N:(x) MG Bs | = [N] {u} u, contribution 


This is the displacement field in terms of nodal variables. 


Truss Element Example 


Figure depicts a tapered elastic bar subjected to an applied tensile load P at one end 
and attached to a fixed support at the other end. The cross-sectional area varies 
linearly from A, at the fixed support at x = 0 to 4/2 at x = L. Calculate the 
displacement of the end of the bar (a) by modeling the bar as a single element 
having cross-sectional area equal to the area of the actual bar at its midpoint along 
the length, (b) using two bar elements of equal length and similarly evaluating the 
area at the midpoint of each, and compare to the exact solution. 


(a) (b) (c) 


Truss Element Example Solution a) 


For a single element, the cross-sectional area is 34,/4 and the element 
“spring constant” and element equation are: 


mL ltd tr] 
L 4 4L L-1 -1]lu, P 


Applying the constraint condition U, — 0, we find for U, as the 
displacement at x = L 


Truss Element Example Solution b) 


Two elements of equal length L/2 with associated nodal displacements. For element 
1, A, = 7Ay8 so 

AIE u TAYE TAYE 

Ex "BL 4L 


ki = 


while for element 2, we have 


SAo A» E SAVE 546 
Aj = — and ka — — = 
8 L 8(L/2) 4L 


Since no load is applied at the center of the bar, the equilibrium equations for the 
system of two elements 1s: 


ki ki 0 U, F, 
—ki ki +k, —k, U, = 0 
0 — ko ka U3 P 


Applying the constraint condition U, = 0 results in 


~ elle | 
—k» kò Us {P 


Truss Element Example Solution b) 


Adding the two equations gives 


c) Exact solution 1.386 


PL 
AoE 


Truss Element Example Solution Comparison 
Deflection 


Exact 
---- One element 
= — Two elements 


u(x) 


0 01 02 03 04 05 06 07 08 09 1.0 
x 
L 


Shape function for interpolated values: u(x) = (1 — x/L)u, + (x/L)u, 


Truss Element Example Solution Comparison 
Stress 


Node 


2 
Exact 
18 r.--- One element 
1.6 - — Two elements z PEERS | 
1.4 
ST S 
Li 12 
| 
` l 
= 0.8 
SIS ： 
0.6 
0.4 
0.2 


0 01 02 03 04 05 06 07 08 09 1.0 


For stress results are much different, discontinuous for FEA and highly dependent on number of nodes 


Beam Element 


The usual assumptions of elementary beam theory are 
applicable here: 


1. The beam is loaded only in the y direction. 


2. Deflections of the beam are small in comparison to the 
characteristic dimensions of the beam. 


3. The material of the beam is linearly elastic, isotropic, and 
homogeneous. The beam is prismatic and the cross section 
has an axis of symmetry in the plane of bending. 


- d 
€ c 1 
Mr A My pi 


Beam Element 


stiffness matrix 


Ji ku Ay : Ky 
"p Ka kn Ez 
Fy Ka Kn kaz 
my ka ka ke 
4x1 4x 4 


£y | [Yr 
ky |] A 
kz | [vs 
ka || 5 
4x1 


The equation shows that the element stiffness matrix for the beam element isa4 X 
4 matrix. This corresponds to the fact that the element exhibits four degrees of 
freedom and that the displacements are not independent (that is, the body is 


continuous and elastic). 


Furthermore, the matrix is symmetric. This is a consequence of the symmetry of the 


forces and moments (equal and opposite to ensure equilibrium). 


Also the matrix is singular and therefore not invertible. That 1s because the 
problem as defined 1s incomplete and does not have a solution: boundary 


conditions are required. 


Beam Element Shape Function and 
Stiffness Matrix 
Shape function: 
v(x) = f(vi, va, 01, 02, x) 
v(x) = (1 — 38% + 2E” vi + L(E — 2E” + E501 + (38? — 2E w 
+ LEE — 10; 


h € ㆍ 
Wit = — 
À L 


And the Stiffness Matrix: 


W OL. =12 GE 

El.| 6L 4L? —6L OLI? 
~ P|- 6. Y 6L 
6L 2L? -—6L AIL? 


Way of Stacking Blocks/Elements 


e Compatibility requirement: ensures that the 
“displacements” at the shared node of adjacent 
elements are equal. 

* Equilibrium requirement: ensures that elemental 
forces and the external forces applied to the 
system nodes are in equilibrium. 

ㆍ Boundary conditions: ensures the system satisfy 
the boundary constraints and so on. 


Cd:stribibkes Load) 


e Belleville (Conical) Spring: 


Limitations 
of Regular 
FEA 


Software 


* Unable to handle p 
geometrically 
nonlinear - large  . 
deformation problems: 


6 
shells, rubber, etc. 


This is a geometrically nonlinear (large deformation) 
problem and iteration method (incremental approach) needs to 
be employed. 


Interpolation Functions for General 
Element Formulation 


In finite element analysis, solution accuracy 1s judged in terms of 
convergence as the element “mesh” 1s refined. 
There are two major methods of mesh refinement. 


In the first, known as h-refinement, mesh refinement refers to the process 
of increasing the number of elements used to model a given domain, 
consequently, reducing individual element size. 


In the second method, p-refinement, element size 1s unchanged but the 
order of the polynomials used as interpolation functions 1s increased. 


The objective of mesh refinement in either method is to obtain sequential 
solutions that exhibit asymptotic convergence to values representing the 
exact solution. 


